Application Final Version

Import packages

library(sf)
library(mapview)
library(dplyr)
library(tidyr)
library(lubridate)
library(evd)  
library(parallel)  
library(FNN)
library(stringr)

Load map

pathshp <- "~/comp9991/data/geo/DOPGrid.shp"
map <- st_read(pathshp, quiet = TRUE)
print(map)

Simple feature collection with 17450 features and 4 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 274189.9 ymin: 2901792 xmax: 532797 ymax: 3324561
Projected CRS: NAD83 / UTM zone 17N
First 10 features:
   OBJECTID  PIXEL Shape_Leng Shape_Area                       geometry
1     14202 131641   7697.982    3699323 POLYGON ((411365.9 3234204,...
2     14203 131642   7697.949    3699292 POLYGON ((413224.1 3234190,...
3     14204 131643   7697.920    3699264 POLYGON ((415082.4 3234176,...
4     14205 131644   7697.892    3699237 POLYGON ((416940.6 3234162,...
5     14206 131645   7697.862    3699208 POLYGON ((418798.9 3234149,...
6     14207 131646   7697.834    3699181 POLYGON ((420657.1 3234136,...
7     14208 131647   7697.806    3699155 POLYGON ((422515.3 3234123,...
8     14209 131648   7697.780    3699130 POLYGON ((424373.5 3234111,...
9     14210 131649   7697.756    3699106 POLYGON ((426231.8 3234098,...
10    14211 131650   7697.728    3699080 POLYGON ((428090 3234087, 4...

mapview(map)

Load data

data <- read.csv("~/comp9991/data/PixelRain15min_1995.csv",header = TRUE)

change_time_resolution <- function(data, res = "default"){
  # data: dataframe
  data <- data %>%
  mutate(DateTime = dmy_hms(DateTime),    
         Date = as.Date(DateTime),
         Hour = floor_date(DateTime, "hour"))
  
  if (res == "default"){
    return(data)
  }
  else if (res == "hourly"){
    data_hourly <- data %>% 
      group_by(Hour) %>% 
      summarise(across(starts_with("P"), \(x) sum(x, na.rm = TRUE)))
    
    return(data_hourly)
  }
  else if (res == "daily"){
    data_daily <- data %>%
      group_by(Date) %>%
      summarise(across(starts_with("P"), \(x) sum(x, na.rm = TRUE)))
    return(data_daily)
  }
}

hourly

data_hourly <- change_time_resolution(data,res = "hourly")
class(data_hourly)

[1] "tbl_df"     "tbl"        "data.frame"

all_pixel_IDs <- names(data_hourly)[-1]
all_pixel_IDs <- as.numeric(stringi::stri_extract_first(all_pixel_IDs, regex = "[0-9]+"))

sub_map <- map %>%
  filter(PIXEL %in% all_pixel_IDs)

mapview(sub_map)

cut a square: 14*14

cutid1 <- seq(11951+2, 11951 + 13+2)
cutid2 <- seq(11874+2, 11874 + 13+2)
cutid3 <- seq(11796+2, 11796 + 13+2)
cutid4 <- seq(11716+2, 11716 + 13+2)
cutid5 <- seq(11635+2, 11635 + 13+2)
cutid6 <- seq(11553+2, 11553 + 13+2)
cutid7 <- seq(11470+2, 11470 + 13+2)
cutid8 <- seq(11386+2, 11386 + 13+2)
cutid9 <- seq(11301+2, 11301 + 13+2)
cutid10 <- seq(11215+2, 11215 + 13+2)
cutid11 <- seq(11129+2, 11129 + 13+2)
cutid12 <- seq(11042+2, 11042 + 13+2)
cutid13 <- seq(10954+2, 10954 + 13+2)
cutid14 <- seq(10866+2, 10866 + 13+2)

cut_id <- c(
  cutid1,
  cutid2,
  cutid3,
  cutid4,
  cutid5, 
  cutid6,
  cutid7,
  cutid8,
  cutid9,
  cutid10,
  cutid11,
  cutid12,
  cutid13,
  cutid14
)

square_map <- sub_map %>%
  filter(OBJECTID %in% cut_id)

mapview(square_map)

data_hourly <- data_hourly[rowSums(data_hourly[, -1]) > 0, ]
sub_pixel <- square_map$PIXEL
idx <- all_pixel_IDs %in% sub_pixel
idx <- c(TRUE, idx)
data_hourly_sub <- data_hourly[, idx]

val <- data_hourly_sub[20, 2:ncol(data_hourly_sub)]

val_long <- val %>%
  pivot_longer(
    cols = starts_with("P"),
    names_to = "PIXEL_COL_NAME",
    values_to = "value"  
  ) %>%
  mutate(
    PIXEL = as.integer(str_remove(PIXEL_COL_NAME, pattern = "P"))
  ) %>%
  select(PIXEL, value)

square_map_value <- square_map %>%
  left_join(val_long, by = "PIXEL")

mapview(square_map_value, zcol="value")

filter extremes

rain_hourly_mat <- as.matrix(data_hourly_sub[, -1])
rain_hourly_uniform <- apply(rain_hourly_mat, 2, function(col) {
  ind <- col > 0
  u <- numeric(length(col))
  u[ind] <- rank(col[ind]) / (sum(ind) + 1)
  u
})

eps <- 1e-10
rain_hourly_uniform[rain_hourly_uniform <= 0] <- eps
rain_hourly_uniform[rain_hourly_uniform >= 1] <- 1 - eps

rain_hourly_pareto <- qgpd(rain_hourly_uniform, xi = 1, beta = 1, mu = 1)

L1_risk_hourly <- rowMeans(rain_hourly_pareto, na.rm = TRUE)
u_L1_hourly   <- quantile(L1_risk_hourly, 0.90, na.rm = TRUE)
extreme_L1_hourly <- rain_hourly_pareto[L1_risk_hourly > u_L1_hourly, ]

idx of timestep

idx_ext_hourly <- which(L1_risk_hourly > u_L1_hourly & !is.na(L1_risk_hourly))
times_hourly <- data_hourly_sub[[1]][idx_ext_hourly]

compute the coord

square_map_m <- st_transform(square_map, crs = 3577)
centroids_m <- st_centroid(square_map_m)

Warning: st_centroid assumes attributes are constant over geometries

coords_m <- st_coordinates(centroids_m)    

nn <- get.knn(coords_m, k = 2)$nn.dist[, 2]
cell_len_m <- median(nn)

coord <- coords_m / cell_len_m

bayesian mixture model

build_Gamma <- function(coord, alpha, theta){
  D <- as.matrix(dist(coord))
  G <- (D / theta)^alpha
  diag(G) <- 0 # \gamma(0) = 0
  return(G)
}

prep_br <- function(Gamma) {
  d <- nrow(Gamma)
  out <- vector("list", d)
  for (j in 1:d) {
    idx <- setdiff(seq_len(d), j)
    Sigj <- ( outer(Gamma[idx, j], rep(1, d-1)) +
                outer(rep(1, d-1), Gamma[idx, j]) -
                Gamma[idx, idx] )
    Sigj <- Sigj + .Machine$double.eps * diag(d-1)
    Lj   <- chol(Sigj)
    out[[j]] <- list(idx = idx, L = Lj, Gamma_col = Gamma[, j])
  }
  out
}

intensity_logskew <- function(x,par,alpha.para=TRUE,log=TRUE){
  sigma = par[[1]]
  if(!is.matrix(x)){x <- matrix(x,nrow=1)}
  n = ncol(x)
  if(n==1) return(1/(x^2))
  omega2 = diag(sigma)
  chol.sigma = chol(sigma)
  inv.sigma = chol2inv(chol.sigma)
  logdet.sigma = sum(log(diag(chol.sigma)))*2
  if(alpha.para){
    alpha = par[[2]]
    delta = c(sigma %*% alpha)/sqrt(c(1+alpha %*% sigma %*% alpha))
  }else{
    delta = par[[2]]
    alpha = c(1 - delta %*% inv.sigma %*% delta)^(-1/2) * c(inv.sigma %*% delta)
  }
  a = log(2) + pnorm(delta,log.p=TRUE)
  q = rowSums(inv.sigma)
  sum.q = sum(q);sum.alpha = sum(alpha)
  q.mat = matrix(q,n,n,byrow=TRUE)
  x.log = log(x)
  x.circ = x.log + matrix(a,nrow=nrow(x),ncol=n,byrow=TRUE)
  beta = (1+sum.alpha^2/sum.q)^(-0.5)
  tau.tilde = apply(x.circ,1,function(x.i)  beta * sum((alpha - sum.alpha*q/sum.q) * (x.i + omega2/2))+ beta*sum.alpha/sum.q)
  A = inv.sigma - q %*% t(q)/sum.q
  val = -(n-3)/2 * log(2) - (n-1)/2*log(pi)-1/2*logdet.sigma - 1/2*log(sum.q) - 1/2 * (sum(q*omega2)-1)/sum.q - 1/8*c(omega2 %*% A %*% omega2) 
  val = val - rowSums(x.log) - 1/2 * apply(x.circ,1,function(x.i) c(x.i %*% A %*% x.i) + sum(x.i * (2*q/sum.q + c(A %*% omega2)))) + pnorm(tau.tilde,log.p=TRUE)
  if(log)
    return(val)
  else
    return(exp(val))    
}

dlog_beta  <- function(x, a, b) dbeta(x, a, b, log = TRUE) # mixture weight prior
dlog_gamma <- function(x, a, b) dgamma(x, shape = a, rate = b, log = TRUE)

component_loglik <- function(X, coord, alpha, theta, anchor = 1){
  Gamma <- build_Gamma(coord = coord, alpha = alpha, theta = theta)
  prep  <- prep_br(Gamma)
  idx   <- prep[[anchor]]$idx
  L     <- prep[[anchor]]$L
  Sigma <- t(L) %*% L
  Xsub  <- as.matrix(X[, idx, drop = FALSE])
  intensity_logskew(
    Xsub, par = list(Sigma, delta = rep(0, ncol(Xsub))),
    alpha.para = FALSE, log = TRUE
  )
}

mix_loglik <- function(X, coord, alpha, th1, th2, anchor = 1) {
  if (!is.finite(th1) || th1 <= 0) return(list(l1 = rep(-Inf, nrow(X)), l2 = rep(-Inf, nrow(X))))
  if (!is.finite(th2) || th2 <= 0) return(list(l1 = rep(-Inf, nrow(X)), l2 = rep(-Inf, nrow(X))))
  l1 <- component_loglik(X, coord, alpha = alpha, theta = th1, anchor = anchor)
  l2 <- component_loglik(X, coord, alpha = alpha, theta = th2, anchor = anchor)
  list(l1 = l1, l2 = l2)
}

tau_from_ll <- function(l1, l2, lambda) {
  a <- log(lambda)     + l1
  b <- log(1 - lambda) + l2
  m <- pmax(a, b)
  exp(a - m) / (exp(a - m) + exp(b - m))
}

mcmc_br <- function(
    data, coord, alpha,
    n_iter   = 4000,
    burn     = 1000,
    thin     = 2,
    # priors
    a_lambda = 1, b_lambda = 1,   # for mixture weight lambda
    a_th2 = 1, b_th2 = 1,         # for theta2
    a_del = 1, b_del = 1,         # for delta
    # proposal sds on log-scale
    prop_sd_log_th2  = 0.10,
    prop_sd_log_del  = 0.10,
    # initials
    init_lambda = 0.5,  # mixture weight
    init_th2    = 15,   # theta2
    init_del    = 20,   # delta
    anchor = 1,
    verbose_every = 100
){
  stopifnot(init_lambda > 0, init_lambda < 1, init_th2 > 0, init_del > 0)
  
  keep_idx <- seq.int(burn + 1, n_iter, by = thin)
  n_keep   <- length(keep_idx)
  
  draws <- data.frame(
    iter     = keep_idx,
    lambda   = NA_real_,   # mixture weight
    theta1   = NA_real_,
    theta2   = NA_real_
  )
  z_draws <- matrix(NA_integer_, nrow = n_keep, ncol = nrow(data))
  tau_bar <- rep(0, nrow(data))
  
  # init
  lambda <- init_lambda
  th2    <- init_th2
  del    <- init_del
  th1    <- th2 + del
  
  ll  <- mix_loglik(X = data, coord = coord, alpha = alpha, th1 = th1, th2 = th2, anchor = anchor)
  tau <- tau_from_ll(ll$l1, ll$l2, lambda = lambda)
  z   <- rbinom(n = length(tau), size = 1, prob = tau)
  
  acc_block <- 0L; n_prop <- 0L
  
  # helper: complete-data log-posterior in η-space (η2=log θ2, ηd=log δ)
  logpost_eta <- function(eta2, etad, z_vec) {
    th2p <- exp(eta2)
    delp <- exp(etad)
    th1p <- th2p + delp
    
    llp <- mix_loglik(X = data, coord = coord, alpha = alpha,
                      th1 = th1p, th2 = th2p, anchor = anchor)
    
    lp_ll    <- sum(z_vec * llp$l1 + (1 - z_vec) * llp$l2)
    lp_prior <- dlog_gamma(th2p, a_th2, b_th2) + dlog_gamma(delp, a_del, b_del)
    lp_jac   <- eta2 + etad
    list(value = lp_ll + lp_prior + lp_jac,
         th1 = th1p, th2 = th2p, del = delp, l1_vec = llp$l1, l2_vec = llp$l2)
  }
  
  keep_cursor <- 1L
  for (it in seq_len(n_iter)) {
    
    ## (1) blocked MH for (θ2, δ) on log-scale
    n_prop <- n_prop + 1L
    eta2  <- log(th2);  etad <- log(del)
    eta2p <- rnorm(1, mean = eta2, sd = prop_sd_log_th2)
    etadp <- rnorm(1, mean = etad, sd = prop_sd_log_del)
    
    # current value via same function
    curr <- logpost_eta(eta2, etad, z)
    prop <- logpost_eta(eta2p, etadp, z)
    
    if (log(runif(1)) < (prop$value - curr$value)) {
      th2 <- prop$th2;  del <- prop$del;  th1 <- prop$th1
      ll$l1 <- prop$l1_vec; ll$l2 <- prop$l2_vec
      acc_block <- acc_block + 1L
    }
    # else keep (th1,th2,del,ll) as-is
    
    ## (2) Gibbs for lambda | z (Beta-conjugate)
    n1 <- sum(z); n0 <- length(z) - n1
    lambda <- rbeta(1, a_lambda + n1, b_lambda + n0)
    
    ## (3) Gibbs for z | x, θ, λ
    tau <- tau_from_ll(ll$l1, ll$l2, lambda)
    z   <- rbinom(length(tau), 1, tau)
    
    ## store thinned draws
    if (it %in% keep_idx) {
      tau_bar <- tau_bar + tau / n_keep
      draws$lambda[keep_cursor] <- lambda
      draws$theta1[keep_cursor] <- th1
      draws$theta2[keep_cursor] <- th2
      z_draws[keep_cursor, ]    <- z
      keep_cursor <- keep_cursor + 1L
    }
    
    if (verbose_every > 0 && it %% verbose_every == 0) {
      cat(sprintf("iter %d | lambda=%.3f θ1=%.3f θ2=%.3f | acc(θ-block)=%.3f\n",
                  it, lambda, th1, th2, acc_block / n_prop))
    }
  }
  
  list(
    draws = draws,
    z_draws = z_draws,
    accept_rate_theta_block = acc_block / n_prop,
    tau_bar = tau_bar,
    hard_labels = ifelse(tau_bar > 0.5, 1, 2),
    post_mean = c(lambda = mean(draws$lambda),
                  theta1 = mean(draws$theta1),
                  theta2 = mean(draws$theta2)),
    post_median = c(lambda = median(draws$lambda),
                    theta1 = median(draws$theta1),
                    theta2 = median(draws$theta2))
  )
}

fit the data

alpha <- 1.5

fit1 <- mcmc_br(
  data = extreme_L1_hourly, coord = coord, alpha = alpha,
  n_iter = 10000, burn = 2000, thin = 2,
  # priors
  a_lambda = 1, b_lambda = 1,   # prior for mixture weight lambda
  a_th2 = 1, b_th2 = 1,         # prior for theta2
  a_del = 1, b_del = 1,         # prior for delta
  # initials
  init_th2  = 2,
  init_del  = 4,
  # proposals (log-scale)
  prop_sd_log_th2  = 0.01,
  prop_sd_log_del  = 0.01,
  anchor = 1,
  verbose_every = 50
)

iter 50 | lambda=0.429 θ1=6.866 θ2=2.087 | acc(θ-block)=0.440
iter 100 | lambda=0.419 θ1=6.927 θ2=2.082 | acc(θ-block)=0.330
iter 150 | lambda=0.421 θ1=6.854 θ2=2.088 | acc(θ-block)=0.287
iter 200 | lambda=0.402 θ1=6.889 θ2=2.100 | acc(θ-block)=0.265
iter 250 | lambda=0.379 θ1=6.896 θ2=2.104 | acc(θ-block)=0.276
iter 300 | lambda=0.393 θ1=6.877 θ2=2.098 | acc(θ-block)=0.297
iter 350 | lambda=0.390 θ1=6.863 θ2=2.081 | acc(θ-block)=0.300
iter 400 | lambda=0.422 θ1=6.822 θ2=2.092 | acc(θ-block)=0.285
iter 450 | lambda=0.363 θ1=6.900 θ2=2.092 | acc(θ-block)=0.280
iter 500 | lambda=0.417 θ1=6.800 θ2=2.094 | acc(θ-block)=0.274
iter 550 | lambda=0.428 θ1=6.818 θ2=2.084 | acc(θ-block)=0.291
iter 600 | lambda=0.376 θ1=6.847 θ2=2.092 | acc(θ-block)=0.288
iter 650 | lambda=0.447 θ1=6.826 θ2=2.092 | acc(θ-block)=0.275
iter 700 | lambda=0.397 θ1=6.877 θ2=2.083 | acc(θ-block)=0.283
iter 750 | lambda=0.388 θ1=6.921 θ2=2.087 | acc(θ-block)=0.285
iter 800 | lambda=0.378 θ1=6.824 θ2=2.088 | acc(θ-block)=0.286
iter 850 | lambda=0.411 θ1=6.806 θ2=2.087 | acc(θ-block)=0.287
iter 900 | lambda=0.426 θ1=6.891 θ2=2.106 | acc(θ-block)=0.289
iter 950 | lambda=0.444 θ1=6.879 θ2=2.102 | acc(θ-block)=0.288
iter 1000 | lambda=0.428 θ1=6.795 θ2=2.092 | acc(θ-block)=0.288
iter 1050 | lambda=0.383 θ1=6.841 θ2=2.084 | acc(θ-block)=0.293
iter 1100 | lambda=0.426 θ1=6.824 θ2=2.083 | acc(θ-block)=0.292
iter 1150 | lambda=0.382 θ1=6.842 θ2=2.098 | acc(θ-block)=0.295
iter 1200 | lambda=0.419 θ1=6.878 θ2=2.094 | acc(θ-block)=0.293
iter 1250 | lambda=0.400 θ1=6.866 θ2=2.111 | acc(θ-block)=0.296
iter 1300 | lambda=0.393 θ1=6.793 θ2=2.089 | acc(θ-block)=0.293
iter 1350 | lambda=0.408 θ1=6.819 θ2=2.099 | acc(θ-block)=0.293
iter 1400 | lambda=0.413 θ1=6.850 θ2=2.090 | acc(θ-block)=0.286
iter 1450 | lambda=0.425 θ1=6.883 θ2=2.094 | acc(θ-block)=0.291
iter 1500 | lambda=0.418 θ1=6.905 θ2=2.086 | acc(θ-block)=0.293
iter 1550 | lambda=0.396 θ1=6.833 θ2=2.091 | acc(θ-block)=0.291
iter 1600 | lambda=0.415 θ1=6.790 θ2=2.091 | acc(θ-block)=0.292
iter 1650 | lambda=0.421 θ1=6.896 θ2=2.081 | acc(θ-block)=0.293
iter 1700 | lambda=0.404 θ1=6.905 θ2=2.093 | acc(θ-block)=0.291
iter 1750 | lambda=0.407 θ1=6.872 θ2=2.084 | acc(θ-block)=0.295
iter 1800 | lambda=0.429 θ1=6.812 θ2=2.103 | acc(θ-block)=0.297
iter 1850 | lambda=0.405 θ1=6.776 θ2=2.095 | acc(θ-block)=0.299
iter 1900 | lambda=0.431 θ1=6.843 θ2=2.088 | acc(θ-block)=0.300
iter 1950 | lambda=0.434 θ1=6.839 θ2=2.095 | acc(θ-block)=0.298
iter 2000 | lambda=0.416 θ1=6.843 θ2=2.085 | acc(θ-block)=0.296
iter 2050 | lambda=0.422 θ1=6.875 θ2=2.093 | acc(θ-block)=0.298
iter 2100 | lambda=0.401 θ1=6.827 θ2=2.089 | acc(θ-block)=0.298
iter 2150 | lambda=0.412 θ1=6.935 θ2=2.082 | acc(θ-block)=0.297
iter 2200 | lambda=0.426 θ1=6.834 θ2=2.096 | acc(θ-block)=0.298
iter 2250 | lambda=0.370 θ1=6.881 θ2=2.090 | acc(θ-block)=0.299
iter 2300 | lambda=0.465 θ1=6.816 θ2=2.096 | acc(θ-block)=0.297
iter 2350 | lambda=0.427 θ1=6.869 θ2=2.088 | acc(θ-block)=0.297
iter 2400 | lambda=0.450 θ1=6.830 θ2=2.095 | acc(θ-block)=0.296
iter 2450 | lambda=0.444 θ1=6.794 θ2=2.086 | acc(θ-block)=0.294
iter 2500 | lambda=0.412 θ1=6.773 θ2=2.082 | acc(θ-block)=0.294
iter 2550 | lambda=0.424 θ1=6.824 θ2=2.107 | acc(θ-block)=0.294
iter 2600 | lambda=0.387 θ1=6.920 θ2=2.109 | acc(θ-block)=0.293
iter 2650 | lambda=0.413 θ1=6.863 θ2=2.099 | acc(θ-block)=0.295
iter 2700 | lambda=0.444 θ1=6.824 θ2=2.091 | acc(θ-block)=0.295
iter 2750 | lambda=0.397 θ1=6.922 θ2=2.106 | acc(θ-block)=0.293
iter 2800 | lambda=0.375 θ1=6.832 θ2=2.100 | acc(θ-block)=0.297
iter 2850 | lambda=0.409 θ1=6.829 θ2=2.091 | acc(θ-block)=0.296
iter 2900 | lambda=0.438 θ1=6.822 θ2=2.081 | acc(θ-block)=0.298
iter 2950 | lambda=0.435 θ1=6.937 θ2=2.100 | acc(θ-block)=0.298
iter 3000 | lambda=0.472 θ1=6.796 θ2=2.085 | acc(θ-block)=0.299
iter 3050 | lambda=0.402 θ1=6.873 θ2=2.094 | acc(θ-block)=0.299
iter 3100 | lambda=0.437 θ1=6.870 θ2=2.094 | acc(θ-block)=0.298
iter 3150 | lambda=0.445 θ1=6.828 θ2=2.093 | acc(θ-block)=0.299
iter 3200 | lambda=0.420 θ1=6.910 θ2=2.088 | acc(θ-block)=0.299
iter 3250 | lambda=0.431 θ1=6.842 θ2=2.097 | acc(θ-block)=0.298
iter 3300 | lambda=0.411 θ1=6.885 θ2=2.088 | acc(θ-block)=0.299
iter 3350 | lambda=0.394 θ1=6.894 θ2=2.098 | acc(θ-block)=0.299
iter 3400 | lambda=0.395 θ1=6.814 θ2=2.107 | acc(θ-block)=0.299
iter 3450 | lambda=0.380 θ1=6.784 θ2=2.092 | acc(θ-block)=0.300
iter 3500 | lambda=0.484 θ1=6.853 θ2=2.087 | acc(θ-block)=0.300
iter 3550 | lambda=0.384 θ1=6.856 θ2=2.095 | acc(θ-block)=0.300
iter 3600 | lambda=0.427 θ1=6.850 θ2=2.084 | acc(θ-block)=0.301
iter 3650 | lambda=0.404 θ1=6.902 θ2=2.097 | acc(θ-block)=0.300
iter 3700 | lambda=0.423 θ1=6.773 θ2=2.089 | acc(θ-block)=0.299
iter 3750 | lambda=0.459 θ1=6.865 θ2=2.110 | acc(θ-block)=0.300
iter 3800 | lambda=0.406 θ1=6.808 θ2=2.088 | acc(θ-block)=0.301
iter 3850 | lambda=0.414 θ1=6.800 θ2=2.097 | acc(θ-block)=0.300
iter 3900 | lambda=0.421 θ1=6.834 θ2=2.097 | acc(θ-block)=0.300
iter 3950 | lambda=0.432 θ1=6.845 θ2=2.104 | acc(θ-block)=0.300
iter 4000 | lambda=0.431 θ1=6.840 θ2=2.095 | acc(θ-block)=0.299
iter 4050 | lambda=0.412 θ1=6.800 θ2=2.090 | acc(θ-block)=0.300
iter 4100 | lambda=0.409 θ1=6.854 θ2=2.096 | acc(θ-block)=0.299
iter 4150 | lambda=0.392 θ1=6.789 θ2=2.090 | acc(θ-block)=0.300
iter 4200 | lambda=0.414 θ1=6.844 θ2=2.097 | acc(θ-block)=0.300
iter 4250 | lambda=0.434 θ1=6.853 θ2=2.077 | acc(θ-block)=0.300
iter 4300 | lambda=0.387 θ1=6.849 θ2=2.095 | acc(θ-block)=0.299
iter 4350 | lambda=0.414 θ1=6.901 θ2=2.097 | acc(θ-block)=0.297
iter 4400 | lambda=0.409 θ1=6.923 θ2=2.086 | acc(θ-block)=0.297
iter 4450 | lambda=0.413 θ1=6.781 θ2=2.091 | acc(θ-block)=0.297
iter 4500 | lambda=0.420 θ1=6.883 θ2=2.095 | acc(θ-block)=0.298
iter 4550 | lambda=0.437 θ1=6.876 θ2=2.091 | acc(θ-block)=0.299
iter 4600 | lambda=0.375 θ1=6.823 θ2=2.099 | acc(θ-block)=0.298
iter 4650 | lambda=0.438 θ1=6.866 θ2=2.091 | acc(θ-block)=0.298
iter 4700 | lambda=0.420 θ1=6.817 θ2=2.094 | acc(θ-block)=0.299
iter 4750 | lambda=0.412 θ1=6.921 θ2=2.094 | acc(θ-block)=0.299
iter 4800 | lambda=0.391 θ1=6.841 θ2=2.093 | acc(θ-block)=0.299
iter 4850 | lambda=0.402 θ1=6.868 θ2=2.100 | acc(θ-block)=0.299
iter 4900 | lambda=0.390 θ1=6.861 θ2=2.100 | acc(θ-block)=0.300
iter 4950 | lambda=0.407 θ1=6.825 θ2=2.090 | acc(θ-block)=0.300
iter 5000 | lambda=0.395 θ1=6.845 θ2=2.096 | acc(θ-block)=0.299
iter 5050 | lambda=0.407 θ1=6.833 θ2=2.090 | acc(θ-block)=0.298
iter 5100 | lambda=0.396 θ1=6.895 θ2=2.091 | acc(θ-block)=0.298
iter 5150 | lambda=0.368 θ1=6.906 θ2=2.109 | acc(θ-block)=0.297
iter 5200 | lambda=0.413 θ1=6.820 θ2=2.086 | acc(θ-block)=0.299
iter 5250 | lambda=0.426 θ1=6.832 θ2=2.099 | acc(θ-block)=0.298
iter 5300 | lambda=0.457 θ1=7.003 θ2=2.097 | acc(θ-block)=0.298
iter 5350 | lambda=0.407 θ1=6.853 θ2=2.087 | acc(θ-block)=0.298
iter 5400 | lambda=0.396 θ1=6.839 θ2=2.095 | acc(θ-block)=0.298
iter 5450 | lambda=0.373 θ1=6.909 θ2=2.092 | acc(θ-block)=0.298
iter 5500 | lambda=0.405 θ1=6.892 θ2=2.080 | acc(θ-block)=0.298
iter 5550 | lambda=0.404 θ1=6.831 θ2=2.089 | acc(θ-block)=0.298
iter 5600 | lambda=0.442 θ1=6.887 θ2=2.088 | acc(θ-block)=0.297
iter 5650 | lambda=0.386 θ1=6.843 θ2=2.078 | acc(θ-block)=0.297
iter 5700 | lambda=0.393 θ1=6.896 θ2=2.096 | acc(θ-block)=0.297
iter 5750 | lambda=0.396 θ1=6.903 θ2=2.101 | acc(θ-block)=0.297
iter 5800 | lambda=0.429 θ1=6.786 θ2=2.094 | acc(θ-block)=0.297
iter 5850 | lambda=0.437 θ1=6.881 θ2=2.083 | acc(θ-block)=0.297
iter 5900 | lambda=0.443 θ1=6.841 θ2=2.098 | acc(θ-block)=0.297
iter 5950 | lambda=0.409 θ1=6.806 θ2=2.096 | acc(θ-block)=0.297
iter 6000 | lambda=0.426 θ1=6.859 θ2=2.097 | acc(θ-block)=0.297
iter 6050 | lambda=0.379 θ1=6.868 θ2=2.092 | acc(θ-block)=0.298
iter 6100 | lambda=0.402 θ1=6.841 θ2=2.103 | acc(θ-block)=0.298
iter 6150 | lambda=0.427 θ1=6.837 θ2=2.093 | acc(θ-block)=0.299
iter 6200 | lambda=0.413 θ1=6.866 θ2=2.069 | acc(θ-block)=0.300
iter 6250 | lambda=0.425 θ1=6.815 θ2=2.083 | acc(θ-block)=0.300
iter 6300 | lambda=0.459 θ1=6.787 θ2=2.097 | acc(θ-block)=0.300
iter 6350 | lambda=0.415 θ1=6.884 θ2=2.107 | acc(θ-block)=0.300
iter 6400 | lambda=0.355 θ1=6.890 θ2=2.096 | acc(θ-block)=0.300
iter 6450 | lambda=0.413 θ1=6.852 θ2=2.079 | acc(θ-block)=0.300
iter 6500 | lambda=0.429 θ1=6.851 θ2=2.105 | acc(θ-block)=0.301
iter 6550 | lambda=0.449 θ1=6.868 θ2=2.075 | acc(θ-block)=0.300
iter 6600 | lambda=0.399 θ1=6.870 θ2=2.098 | acc(θ-block)=0.300
iter 6650 | lambda=0.432 θ1=6.878 θ2=2.107 | acc(θ-block)=0.300
iter 6700 | lambda=0.423 θ1=6.850 θ2=2.099 | acc(θ-block)=0.300
iter 6750 | lambda=0.437 θ1=6.833 θ2=2.099 | acc(θ-block)=0.299
iter 6800 | lambda=0.386 θ1=6.915 θ2=2.101 | acc(θ-block)=0.300
iter 6850 | lambda=0.443 θ1=6.844 θ2=2.093 | acc(θ-block)=0.300
iter 6900 | lambda=0.398 θ1=6.832 θ2=2.103 | acc(θ-block)=0.299
iter 6950 | lambda=0.385 θ1=6.863 θ2=2.100 | acc(θ-block)=0.299
iter 7000 | lambda=0.353 θ1=6.860 θ2=2.104 | acc(θ-block)=0.299
iter 7050 | lambda=0.455 θ1=6.841 θ2=2.089 | acc(θ-block)=0.299
iter 7100 | lambda=0.396 θ1=6.833 θ2=2.095 | acc(θ-block)=0.298
iter 7150 | lambda=0.457 θ1=6.860 θ2=2.080 | acc(θ-block)=0.299
iter 7200 | lambda=0.412 θ1=6.890 θ2=2.101 | acc(θ-block)=0.299
iter 7250 | lambda=0.403 θ1=6.835 θ2=2.089 | acc(θ-block)=0.299
iter 7300 | lambda=0.458 θ1=6.819 θ2=2.083 | acc(θ-block)=0.299
iter 7350 | lambda=0.385 θ1=6.919 θ2=2.083 | acc(θ-block)=0.300
iter 7400 | lambda=0.426 θ1=6.829 θ2=2.109 | acc(θ-block)=0.299
iter 7450 | lambda=0.454 θ1=6.874 θ2=2.100 | acc(θ-block)=0.299
iter 7500 | lambda=0.460 θ1=6.803 θ2=2.097 | acc(θ-block)=0.299
iter 7550 | lambda=0.413 θ1=6.808 θ2=2.104 | acc(θ-block)=0.298
iter 7600 | lambda=0.453 θ1=6.820 θ2=2.092 | acc(θ-block)=0.298
iter 7650 | lambda=0.386 θ1=6.839 θ2=2.098 | acc(θ-block)=0.298
iter 7700 | lambda=0.418 θ1=6.844 θ2=2.098 | acc(θ-block)=0.299
iter 7750 | lambda=0.376 θ1=6.866 θ2=2.082 | acc(θ-block)=0.300
iter 7800 | lambda=0.437 θ1=6.842 θ2=2.083 | acc(θ-block)=0.300
iter 7850 | lambda=0.366 θ1=6.843 θ2=2.085 | acc(θ-block)=0.299
iter 7900 | lambda=0.479 θ1=6.838 θ2=2.094 | acc(θ-block)=0.299
iter 7950 | lambda=0.417 θ1=6.795 θ2=2.093 | acc(θ-block)=0.300
iter 8000 | lambda=0.437 θ1=6.837 θ2=2.081 | acc(θ-block)=0.300
iter 8050 | lambda=0.422 θ1=6.945 θ2=2.096 | acc(θ-block)=0.300
iter 8100 | lambda=0.424 θ1=6.814 θ2=2.091 | acc(θ-block)=0.299
iter 8150 | lambda=0.451 θ1=6.873 θ2=2.100 | acc(θ-block)=0.299
iter 8200 | lambda=0.439 θ1=6.843 θ2=2.085 | acc(θ-block)=0.299
iter 8250 | lambda=0.427 θ1=6.859 θ2=2.099 | acc(θ-block)=0.299
iter 8300 | lambda=0.430 θ1=6.893 θ2=2.094 | acc(θ-block)=0.298
iter 8350 | lambda=0.424 θ1=6.869 θ2=2.102 | acc(θ-block)=0.298
iter 8400 | lambda=0.390 θ1=6.865 θ2=2.085 | acc(θ-block)=0.298
iter 8450 | lambda=0.393 θ1=6.830 θ2=2.092 | acc(θ-block)=0.299
iter 8500 | lambda=0.430 θ1=6.809 θ2=2.077 | acc(θ-block)=0.299
iter 8550 | lambda=0.416 θ1=6.837 θ2=2.100 | acc(θ-block)=0.299
iter 8600 | lambda=0.402 θ1=6.826 θ2=2.085 | acc(θ-block)=0.300
iter 8650 | lambda=0.465 θ1=6.829 θ2=2.091 | acc(θ-block)=0.300
iter 8700 | lambda=0.415 θ1=6.816 θ2=2.096 | acc(θ-block)=0.300
iter 8750 | lambda=0.392 θ1=6.839 θ2=2.081 | acc(θ-block)=0.300
iter 8800 | lambda=0.433 θ1=6.839 θ2=2.099 | acc(θ-block)=0.299
iter 8850 | lambda=0.426 θ1=6.795 θ2=2.091 | acc(θ-block)=0.298
iter 8900 | lambda=0.430 θ1=6.912 θ2=2.083 | acc(θ-block)=0.298
iter 8950 | lambda=0.391 θ1=6.844 θ2=2.099 | acc(θ-block)=0.298
iter 9000 | lambda=0.433 θ1=6.830 θ2=2.102 | acc(θ-block)=0.298
iter 9050 | lambda=0.420 θ1=6.875 θ2=2.101 | acc(θ-block)=0.298
iter 9100 | lambda=0.401 θ1=6.876 θ2=2.093 | acc(θ-block)=0.297
iter 9150 | lambda=0.402 θ1=6.886 θ2=2.107 | acc(θ-block)=0.297
iter 9200 | lambda=0.428 θ1=6.846 θ2=2.092 | acc(θ-block)=0.296
iter 9250 | lambda=0.401 θ1=6.897 θ2=2.094 | acc(θ-block)=0.295
iter 9300 | lambda=0.413 θ1=6.791 θ2=2.093 | acc(θ-block)=0.295
iter 9350 | lambda=0.377 θ1=6.864 θ2=2.098 | acc(θ-block)=0.295
iter 9400 | lambda=0.418 θ1=6.830 θ2=2.100 | acc(θ-block)=0.296
iter 9450 | lambda=0.397 θ1=6.901 θ2=2.081 | acc(θ-block)=0.296
iter 9500 | lambda=0.363 θ1=6.867 θ2=2.092 | acc(θ-block)=0.295
iter 9550 | lambda=0.384 θ1=6.833 θ2=2.092 | acc(θ-block)=0.295
iter 9600 | lambda=0.405 θ1=6.873 θ2=2.087 | acc(θ-block)=0.295
iter 9650 | lambda=0.390 θ1=6.851 θ2=2.099 | acc(θ-block)=0.296
iter 9700 | lambda=0.423 θ1=6.871 θ2=2.081 | acc(θ-block)=0.296
iter 9750 | lambda=0.448 θ1=6.902 θ2=2.092 | acc(θ-block)=0.296
iter 9800 | lambda=0.417 θ1=6.867 θ2=2.097 | acc(θ-block)=0.296
iter 9850 | lambda=0.419 θ1=6.903 θ2=2.093 | acc(θ-block)=0.295
iter 9900 | lambda=0.435 θ1=6.832 θ2=2.102 | acc(θ-block)=0.295
iter 9950 | lambda=0.488 θ1=6.849 θ2=2.091 | acc(θ-block)=0.295
iter 10000 | lambda=0.423 θ1=6.845 θ2=2.100 | acc(θ-block)=0.295

save(fit1, file = "application_fit_hourly_final.RData")

df <- fit1$draws
g1 <- ggplot(df, aes(iter, theta1)) + geom_line() +
  labs(title = expression(trace~theta[1]))
g2 <- ggplot(df, aes(iter, theta2)) + geom_line() +
  labs(title = expression(trace~theta[2]))
g3 <- ggplot(df, aes(iter, lambda)) + geom_line() +
  labs(title = expression(trace~lambda))
gd1 <- ggplot(df, aes(theta1)) + geom_density() +
  labs(title = expression(density~theta[1]))
gd2 <- ggplot(df, aes(theta2)) + geom_density() +
  labs(title = expression(density~theta[2]))
gd3 <- ggplot(df, aes(lambda)) + geom_density() +
  labs(title = expression(density~lambda))

gridExtra::grid.arrange(g1,g2,g3, gd1,gd2,gd3, ncol = 3)

fit1$post_mean

   lambda    theta1    theta2 
0.4186907 6.8560784 2.0937366 

sd(fit1$draws$theta1)

[1] 0.03948329

sd(fit1$draws$theta2)

[1] 0.008061415

plot_extreme_on_map <- function(i, extreme_mat = extreme_L1_hourly,
                                map_sf = square_map, data_wide = data_hourly_sub) {
  stopifnot(is.matrix(extreme_mat),
            i >= 1, i <= nrow(extreme_mat))
  
  # column order of pixels in the wide daily data used to build rain_mat/extreme_L1
  pix_order <- as.integer(sub("P", "", colnames(data_wide)[-1]))
  
  # values for the i-th extreme event
  vals_i <- as.numeric(extreme_mat[i, ])
  df_i <- data.frame(PIXEL = pix_order, val_i = vals_i)
  
  # join to geometry
  joined <- map_sf %>%
    left_join(df_i, by = "PIXEL")
  
  # quick stretch for nicer legend breaks
  brks <- quantile(df_i$val_i, probs = seq(0, 1, length.out = 6), na.rm = TRUE)
  brks_fixed <- sort(unique(brks))
  # plot
  mapview(
    joined,
    zcol       = "val_i",
    layer.name = paste0("Extreme L1 snapshot (row ", i, ")"),
    at         = brks_fixed
  )
}

plot_extreme_on_map(1)

plot_extreme_on_map(5)

\(\chi^2(S_A, S_B)\) test

chi_q <- function(X, q = 0.98) {
  n <- ncol(X)
  m <- nrow(X)
  
  # thresholds at site level
  u <- apply(X, 2, quantile, probs = q, na.rm = TRUE)
  
  # initialize matrix
  chi_hat <- matrix(NA, n, n)
  
  for (i in 1:(n-1)) {
    for (j in (i+1):n) {
      Ai <- X[, i] > u[i]
      Aj <- X[, j] > u[j]
      chi_hat[i, j] <- mean(Ai & Aj) / mean(Ai)
      chi_hat[j, i] <- chi_hat[i, j]
    }
  }
  diag(chi_hat) <- 1
  chi_hat
}

dist_matrix <- as.matrix(dist(coord))

stopifnot(length(fit1$hard_labels) == nrow(extreme_L1_hourly))
labs <- fit1$hard_labels

X1 <- extreme_L1_hourly[labs == 1, , drop = FALSE]
X2 <- extreme_L1_hourly[labs == 2, , drop = FALSE]

# 4) compute χ̂_q for both comps and two q levels
chi_95_c1 <- chi_q(X1, q = 0.95)
chi_95_c2 <- chi_q(X2, q = 0.95)
chi_90_c1 <- chi_q(X1, q = 0.90)
chi_90_c2 <- chi_q(X2, q = 0.90)

upper_to_df <- function(chi_mat, dist_mat, q_label, comp_label) {
  idx <- which(upper.tri(chi_mat), arr.ind = TRUE)
  tibble(
    dist = dist_mat[idx],
    chi  = chi_mat[idx],
    q    = q_label,
    comp = comp_label
  ) %>% filter(is.finite(chi))
}

df_all <- bind_rows(
  upper_to_df(chi_95_c1, dist_matrix, "q=0.95", "Comp1"),
  upper_to_df(chi_95_c2, dist_matrix, "q=0.95", "Comp2"),
  upper_to_df(chi_90_c1, dist_matrix, "q=0.90", "Comp1"),
  upper_to_df(chi_90_c2, dist_matrix, "q=0.90", "Comp2"),

)

# 6) plot (mirrors your original style)
ggplot(df_all, aes(x = dist, y = chi, color = comp)) +
  geom_point(alpha = 0.4, size = 0.6) +
  facet_wrap(~ q, scales = "free_y") +
  theme_bw() +
  labs(
    x = expression(h(s[A], s[B])),
    y = expression(hat(chi)[q](s[A], s[B])),
    title = "Empirical " %+% expression(hat(chi)[q]) %+% " vs Distance by MCMC Hard-Label Component"
  ) +
  scale_color_manual(values = c("Comp1" = "blue", "Comp2" = "red"))

for (i in seq_len(ncol(fit1$z_draws))) {
  plot(
    fit1$z_draws[, i],
    type = "l",
    xlab = "iterations",
    ylab = "z draw",
    main = paste("sample", i) # <--- Modified line
  )
}

plot(times_hourly, fit1$tau_bar, type = "l",
     xlab = "Time", ylab = "Posterior mean of z",
     main = "Posterior mean with season")